function [irf,irfse] = estimate_lp_level(depvar,shock)

irf = NaN(16,1);
irfse = NaN(16,1);

shock = shock./std(shock(~isnan(shock)));

for h=0:16
    y = lagmatrix(depvar,-h);
    X = [shock lagmatrix(shock,1) lagmatrix(depvar,1) ones(size(shock))];
    [beta,Var] = nwest(y,X,h+1);
    irf(h+1)=beta(1);
    irfse(h+1)=sqrt(Var(1,1));
end

end
